Capillary condensation in a square geometry with surface fields 
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We study the influence of wetting on the capillary condensation for a simple fluid in a square 
geometry with surface fields, where the reference system is an infinitely long slit. The Corner 
Transfer Matrix Renormalization Group method has been extended to study a two-dimensional 
Ising model confined in an L x L geometry with equal surface fields. Our results have confirmed 
that in both geometries the coexistence line shift is governed by the same scaling powers, but their 
prefactors are different. 
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Porous materials are solids consisting of interconnected 
network of pores. In recent years, microporous and 
mesoporous materials have been a focus of nanoscience 
and nanotechnology as their properties differ significantly 
from the same bulk materials^i^. Both the pore size and 
its shape, as well as the chemical nature of its surface, 
i.e., whether hydrophilic or hydrophobic, determine the 
properties of porous materials^. In a mesopore, the cu- 
mulative effect of the walls becomes important. After the 
formation of adsorbate layers with 2-3 molecular thick- 
ness on walls, further adsorption induces attractions be- 
tween adsorbate molecules, leading to a sudden conden- 
sation of liquid-like adsorbate molecules inside the pores. 
This effect is analogous to the capillary condensation 
phenomena^. 

The Ising square of a finite size L with the field hi 
acting on all four surfaces (boundaries) can be used as 
an idealized representation of a simple fluid in a pore 
or between finely divided colloidal particles. The choice 
of the model entails the assumption that the intermolec- 
ular forces are short-ranged in character, there are no 
dispersion-like forces. In fluid experiments the equiva- 
lent of the bulk magnetic field H describes the deviation 
of the critical chemical potential i? ^ /i — /iq that is de- 
termined by the density of the fluid in the reservoir. The 
phenomenon equivalent to the capillary condensation can 
be studied in magnetic systems^.. Generally, in the bulk, 
the phase coexistence occurs for temperatures T < 
and for vanishing bulk magnetic field i?. In a slit L x oo 
with identical surface fields at the boundaries, the com- 
bined effect of surface fields and confinement shifts the 
phase coexistence to a non-zero value of the bulk mag- 
netic field H = Hcoe{L), which for large L scales as 



out spins is identified with the wall. The above equation 
is known in the literature as the Kelvin equatiorti. 



Hcoe{L) ~ 



(To cos 9 1 



mb 



(1) 



where a^, nib, and 9 are the surface tension of the free 
up-spin/down-spin interface, bulk spontaneous magne- 
tization, and contact angle given by Young's equation, 
respectively^. In our case the liquid phase is represented 
by up-spins, gas phase by down-spins, whereas area with- 
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FIG. 1. The critical wetting lines for semi-infinite systems 
at r = 2.1: the corner wetting (hcw(T)) on the two planar 
surfaces forming a straight angle and wetting {hiu{T)) on the 
planar surface. With respect to the planar wetting hi =0.1 
and 0.275 values correspond to the dry regime, whereas the 
hi = 0.371 and 0.5 value to the wet one. 

Wetting occurs in systems close to the line of phase 
coexistence when one phase may adsorb preferentially at 
a solid substrate. Typically, we model the substrate using 
a planar surface where the critical line {T) , presented 
in Fig. ([T]), is known exactlyi^ 



exp{2J/kBT) (cosh(2J/fcBr) 



cosh(2/i,„ J/fcsT)) 
= sinh(2 J/ZcsT). (2) 



For this semi-infinite system if we approach the co- 
existence line from the gas side along a given isotherm 
(Tc > T > T^) the amount of liquid adsorbed on the 
surface I diverges I ^ H^^= . Moreover, the phase tran- 
sition, called the complete wetting, is characterized by 
the presence of the singular part of the excess surface 
free-energy fsing = H"^'"" ■ For the two-dimensional 
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Ising model, the values of the critical exponents are 
2-af = 2/3 and /3f = 1/3. 

Consequently, when we consider the (pseudo)two- 
dimensional Ising system in a slit geometry and the in- 
teraction between the liquid and walls is strong (above 
the hu}{T) line), the Kelvin equation fails and some cor- 
rections are necessary. The reason is that, although the 
macroscopically thick layer forms only for a semi-infinite 
system, a noticeable liquid layer intervenes between a gas 
and the wall for a finite-size system as well. Therefore, 
Foster pointed out^ that if adsorbed layers were formed 
prior to condensation, the slit width L in Eq. ([T|) should 
be corrected by the layer thickness I. Derjaguin showed^ 
that if solid-fluid forces decayed exponentially or had a 
finite range, the pore width L could be replaced by L — 21. 

Next, Evans et al.'* showed that the effects of wetting 
layers were of quantitative rather than of qualitative im- 
portance for capillary condensation. Albano et al.''^^ and 
Parry and Evana^^ analyzed the next-order correction 
term to the Kelvin equation for the semi-infinite system. 
Both studies, using scaling and thermodynamics argu- 
ments, concluded that for temperatures below the wet- 
ting temperature (the dry regime) the leading correc- 
tion to scaling term was of type Above the wetting 
(the wet regime) the correction is expected to be non- 
analytic due to a singularity of the surface free energy. 
For the two-dimensional Ising mode l^°'^^ the predicted 
correction term is proportional to L"^/^. 

In the subsequent numerical investigation, the density- 
matrix renormalization techniques were employedi^*^. 
For a large range of surface fields and temperature higher- 
order corrections were not compatible with L~^/^, but 
they were of type L~'^/^. It has been shown that this 
apparent disagreement was due to the fact that even for 
the large sizes considered (L ^ 150) the wetting layer has 
a limited thickness, so that the singular part of the sur- 
face free energy that determines the correction-to-scaling 
behavior is dominated by the contacts with the walls. 

In real systems, where properties of both pure fiuids 
and fiuid mixtures confined to nanoporous and microp- 
orous materials are under consideration^^, the square- like 
geometry is more common than the slit one. In this case, 
the geometry of the system significantly affects the course 
of wetting phenomena because close to a corner the im- 
pact of the individual walls is strong, which should lead to 
more intensive formation of the wetting layer. Therefore, 
the corner wetting transition should also be taken into 
accounti^iil. As two sides of the square form the straight 
angle, the corresponding (L — )■ oo) corner wetting line 
hcw{T), presented in Fig. ([T]), is known exactly^ii^ 



cosh(2 J/fcsT) - exp(-2J/fcBr) sinh^(2 J/fcsT) 

^cosh{2 J hen, I ksT). (3) 

In order to model the influence of wetting phenomena 
on the capillary condensation in a square geometry we 
consider a square Ising ferromagnet subject to identical 



boundary fields with the following Hamiltonian: 
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with J > and Stj = ±1. 



(4) 

The first sum is taken over the 
nearest neighbours, while the second sum is performed 
on spins at the surface only. The surface field hi corre- 
sponding to direct short-range interactions between the 
walls and spins is related to the preferential adsorption 
on the surface for one of the two phases. The uniform 
bulk magnetic field H acts over all spins. 




FIG. 2. Pseudo-coexistence lines for the square (triangles) 
and slit (circles) systems for L = 500 and hi = 0.8. L is 
measured in units of lattice constant, whereas hi is in J units. 

The origin of applied numerical method, called 
the Corner Transfer Matrix Renormalization Group 
(CTMRG), came from Baxter^a. Next Nishino and 
Okunishi^^ combined his corner transfer matrix method 
with the ideas from the density matrix renormalization 
group method (DMRG) approach. The last technique 
was developed by White^^^^"^ for the study of ground state 
properties of quantum spin chains and next extended 
by Nishino to two-dimensional classical systems in a slit 
geometry? 

The general idea is to find a representation of the con- 
figurational space in a restricted space which is much 
smaller than the original one m <^ 2^ . This truncation 
is done through the construction of a reduced density 
matrix whose eigenstates provide the optimal basis set 
TO. Of course, the larger to, the better accuracy, so in 
the present case we keep this parameter up to m = 400. 

Although in the original CTMRG algorithm the full 
transfer matrix is never constructed, we have modified it 
to determine the two eigenvectors related to the largest 
eigenvalues. Because each of these vectors dominates on 
the opposite side of the coexistence line, using both vec- 
tors for the construction of the density matrix guaran- 
tees that the Hamiltonian is properly projected on the 
subspace of most probable states. Our results have not 
shown any ambiguities of the calculated free energy and 
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by increasing the number of states kept m, our results 
(the free energy) converged. To have a point of refer- 
ence, we compared the resuhs for the square geometry 
with the resuhs for the sht geometry, where the DMRG 
technique was applied^ 

Due to the finiteness of L and to the nonvanishing 
surface field hi the (pseudo)coexistence lines (for both 
geometries) are shifted with respect to the bulk coexis- 
tence line {H = 0). As one can see in Fig. ([2|) this effect is 
much stronger for the square-like system. This is easy to 
understand, if one remembers that the influence of walls 
on the individual spins is here reinforced with respect to 
the slit geometry. 

A schematic draw of the fluid shape adopted in a 
square geometry below [hi = 0.1) and above [hi = 0.6) 
wetting is demonstrated in Fig. ([3]). Both upper magne- 
tization profiles correspond to the phase when a liquid 
fills the square-like pore. For a weak field the magneti- 
zation is lower nearby walls, whereas for a strong field 
the magnetization is essentially uniformly high through- 
out the square. The bottom magnetizations illustrate 
the situation on the other side of the coexistence line, 
when the gas phase occupies the middle of the square. 
While the weak-field magnetization profile only slightly 
increases at the square edges, the strong-field magneti- 
zation profile increases considerably, which corresponds 
to creation of the liquid layer. Note that the surface im- 
pact is always enhanced in the square corners, although 
limited to a relatively small area. 
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FIG. 3. Two-dimensional magnetization profiles {si,j) on the 
square geometry 31 x 31 for the two values of the surface 
field hi =0.1 (left graphs) and hi = 0.6 (right graphs) cal- 
culated at fixed temperature T — 2.0. The profiles change 
their shapes and spin polarizations while varying the bulk field 
value H below and above the wetting temperature r,„(/ii — 
0.1) = 2.257 and T^{hi = 0.6) = 1.814, respectively. 

The real phase transition (the complete wetting) oc- 
curs only in a semi-infinite system. Since we deal with 
a finite-size system, the thickness of the wetting layer is 
limited, and a sharp liquid-gas interface is not observed. 



Furthermore, because for a weak surface field (below the 
wetting line) there remains a thin liquid layer between 
a gas and the wall, the gas- wall interface formally never 
occurs. For the same reason, below the wetting line the 
isolated droplets of liquid cannot be observed on the wall, 
so the contact angle cannot be drawn. 

As a useful tool for the analysis of the higher terms of 
the Kelvin equation, we introduce the logarithmic deriva- 
tive which acts as an effective dominant exponent 



a{L) 
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ln{L + AL) - ln(L) 



(5) 



When the following expansion of the Kelvin equation is 
assumed 
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combining both formulae, we obtain the first-order ex- 
pansion for the effective exponent 



a{L) — a + {'^ — a) 



B 1 



(7) 



The Kelvin equation is expected to be valid to the first 
order for all T < Tc, but since only a limited size L is 
available for the numerical computation, it is preferable 
to consider only temperatures not too close to the bulk 
critical temperature, where the scaling of the capillary 
critical point Hcrit{L) ~ L^^/^ is present^. Therefore, 
our results refer to temperature T = 2.1, but we have 
checked that they do not qualitatively change for other 
temperatures. 
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FIG. 4. Plots of the local exponent for both geometries at 
T = 2.1. Lines correspond to the square geometry, whereas 
the symbols are related to the slit one. 

Figure ^ shows that for a wide range of the surface 
fields, when the system grows, the value Oi{L) goes to 
a = 1. It confirms that to the accuracy of the first- 
order expansion, the shift of the phase coexistence and 
the system size are reciprocal in both geometries. 
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For the wet phase {9 = 0) in the sht geometry, the 
exact formula for the coefficient is known (see Eq.([T])) 
giving A{T — 2.1) = 0.335 which very weh agrees 
with our estimated value 0.345. Generally, our numer- 
ical fittings for curves Hcoe{L) according to Eq.® show 
that the ratio of the main prefactors for the slit and 
square geometries is Asm / Asquare ~ 1/2 for each value 
of the surface field. For example, below the wetting 
Asm{hi = 0.1) = 0.121 and Asquare{hi = 0.1) = 0.242, 
whereas above the wetting Asiit{hi — 0.5) — 0.345 and 
Asquare{hi = 0.5) = 0.669. This fact can be explained in 
the following way: when the size of the square increases, 
the impact of the corners remains more or less the same, 
but the influence of the sides increases. Thus, when L 
goes to infinity, the system begins to resemble a set of 
two slits perpendicular to each other. 
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FIG. 5. The effective exponents oi(V) and their fitting curves 
denoted by the thick lines for two sample values of the surface 
field: h\ — 0.1 for the dry system and hi = 0.5 for the wet 
system (see Fig. Q). 

As one can see in Figs. H] and [SJ the shape of the curves 
indicates that we are dealing with various expressions of a 
higher order. First of all, the B coefficient that is always 
positive in the slit geometry, is negative in the dry regime 
for the square geometry. This is manifested by the a{L) 
function minima for small 1/L, where the leading cor- 
rection B/L'' dominates. Obviously, the next correction 
has to be positive. Moreover, in the wet regime, where 
the B coefficients are positive, as one can see in Fig. (jH), 
the curves corresponding to the same value of the surface 
field start to overlap when L becomes enough large. 

In the square geometry, as far as the leading correction 
terms are concerned, the precise analysis of our numerical 
results shows, in both the dry and wet regimes, they are 
the same as in the slit geometry (see Fig. [SJ . 

It is worth adding that the value of the surface field 
hi =0.1 {hi — 0.5) was chosen to be in the dry (wet) 



regime with respect to both wetting curves (see Fig. [Ij . 
Although the presented curves correspond to infinite sys- 
tems, we are dealing with systems large enough to guar- 
antee that both points are localised on the appropriate 
side of the wetting curves. 

To sum up, the Corner Transfer Matrix Renormal- 
ization Group method has been extended to study the 
equilibrium statistical mechanics of simple fluids con- 
fined in the square geometry providing some input for 
understanding experiments in porous materials. As the 
method is not perturbative, it can be applied to arbi- 
trary values of the model parameters and yields accurate 
results, provided that a convergence of the free energy 
is reached. Accuracy of the CTMRG method is com- 
pletely controlled by varying the number of the states 
kept. Moreover, the method does not suffer from getting 
stuck in a local minimum of the free energy instead of 
the true global free energy minimum what is a typical 
problem in Monte Carlo simulations. There are neither 
metastabilities nor hysteresis effects provided that m is 
set to be sufficiently large which is our case concluded by 
the converged free energy. 

Our results confirm that in a more realistic (square) 
geometry the coexistence line shift is inversely propor- 
tional to the system size L and the main pre-factor is 
two times larger than for the case of an infinitely long 
slit of width L. Moreover, for the square geometry the 
leading correction term to the Kelvin equation becomes 
negative in the dry regime. 

We have also found that similarly to the slit geometry 
the leading corrections to scaling in the square geometry 
are of type 1/L^ in the dry regime and of type 1 /L^/^ in 
the wet regime. The last correction is again non-analytic 
due to a singularity of the surface free energy. 

Curve fitting was carried out for the values of the sur- 
face fields that correspond to the dry/ wet regime with re- 
spect to both wetting lines defined by Eqs.[2]and[31 The 
area between the wetting lines was not analysed in de- 
tail, although the strongest competition can be expected 
between the two types of wetting. This requires an ex- 
amination of larger systems, which implies the need for 
a more precise calculation. In our next studies, we aim 
to increase the precision of the numerical results by ad- 
ditional improvement of the CTMRG algorithm. 
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